/**
 * Algorithms for AdaptiveArray.
 *
 * \file aa_algo.hpp
 * \author tomby
 * \date 04-01-2009
 **/

/*  Copyright (C) 2009  Tomas 'tomby' Bily

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>. */

#ifndef __DBC_AA_ALGO__
#define __DBC_AA_ALGO__

#include "shapes.hpp"

namespace dbc {

  template<typename aa_type>
  bool
  coarse_if_possible (typename aa_type::node_iterator& it, int rel_level = 1)
  {
    if (!it.level ())
      return false;

    typename aa_type::level_iterator lit (it.level_iter (rel_level));

    for (; !lit.end_p (); ++lit)
      {
	typename aa_type::NumericArrayType suba (*lit);
	
	if (!aa_type::AdaptiveConstrain::coarse_p (suba))
	  return false;
      }

    it.coarse ();

    return true;
  }

  template<typename aa_type, typename border_iterator>
  void 
  aa_fill (aa_type* array,
	   border_iterator& begin,
	   border_iterator& end,
	   typename aa_type::IndexType& interior_point, 
	   typename aa_type::NumericType val)
  { 
    typedef std::map<typename aa_type::IndexType, bool, axlex_less_index<int, aa_type::Dim> > MapType;
    typedef std::map<typename aa_type::NumericArrayType*, bool> CacheType;
    typedef blitz::TinyVector<int, aa_type::Dim> IndexType;
    typedef typename aa_type::cell_index PairType;

    MapType border_nodes;

    for (border_iterator i = begin; i != end; ++i)
      {
	typename aa_type::node_iterator tmp_it (array->at (*i));
	tmp_it.refine (array->refine_level ());
	if (border_nodes.find (tmp_it.position ()) == border_nodes.end ())
	  border_nodes [tmp_it.position ()] = true;
	(*array) (*i) = val;
      }

    typename aa_type::NumericArrayType tmp (IndexType (1 << array->refine_level ()));
    tmp = val;
    typename aa_type::Node tmp_node (array->refine_level (), tmp);

    std::deque<PairType> it_queue;
    blitz::RectDomain<aa_type::Dim> dom (-1, 1);

    {
      typename aa_type::node_iterator it = array->at (interior_point);
      it.refine (array->refine_level ());
      it_queue.push_back (PairType (it.position (), interior_point));
    }
    
    int tmp_cnt = 0;

    typename aa_type::IndexType s (1 << array->refine_level ());
    CacheType cache;

    while (!it_queue.empty ())
      {
	tmp_cnt ++;

	PairType pair (it_queue.back ());
	it_queue.pop_back ();

	if ((*array) (pair.node(), pair.val()) == val)
	  continue;

	typename aa_type::node_iterator nit = array->at (pair.node(), false);
	typename aa_type::NumericArrayType* data = *nit;

	if (data)
	  {
	    if (cache.find (data) == cache.end ())
	      {
		cache [data] = blitz::all ((*data) != val);
		printf ("cached\n");
	      }
	  }

	if (!data || (data && cache [data]))
	  {
	    nit.coarse (array->refine_level ());
	    nit.data () = val;

	    HyperCubeBorderIterator<aa_type::Dim> hcit (s);

	    nit.node () = &tmp_node;
	    
	    for (; !hcit.end_p (); ++hcit)
	      {
		AANeighbourIterator<aa_type> nit1 (nit, *hcit, dom);

		for (; !nit1.end_p(); ++nit1)
		  if (*nit1 != val)
		    {
		      PairType tmp_p (nit1.global_position (), nit1.local_position ());
		      it_queue.push_back (tmp_p);
		    }
	      }

	    nit.node () = NULL;
	  }
	else
	  {
	    (*array) (pair.node(), pair.val()) = val;
	    
	    AAVirtNeighbourIterator<aa_type> it (nit, &tmp_node, pair.val(), dom);
	    
	    for (; !it.end_p (); ++it)
	      if (*it != val)
		{
		  PairType tmp_p (it.global_position (), it.local_position ());
		  it_queue.push_back (tmp_p);
		}
	  }
	
      }

    printf ("cnt %d\n", tmp_cnt);
#if 1
    for (typename MapType::iterator i = border_nodes.begin (); i != border_nodes.end (); ++i)
      {
	typename aa_type::node_iterator tmp_ni (array->at ((*i).first, false));
	
	while (coarse_if_possible<aa_type> (tmp_ni))
	  ;
      }
#endif
  }

  template<int dim>
  inline blitz::TinyVector<int, dim>
  aa_convert_to_level_coord (blitz::TinyVector<int, dim>& pos, int dest_level, int src_level, bool lb)
  {
    int rel_level = dest_level - src_level;

    if (!rel_level)
      return pos;

    if (rel_level > 0)
      {
	if (lb)
	  return pos << rel_level;
	else
	  return (pos << rel_level) + (1 << rel_level) - 1;
      }

    if (lb)
      return pos >> (-rel_level);
    
    blitz::TinyVector<int, dim> res;

    for (int i = 0; i < dim; i ++)
      if (pos [i] & ((1 << -rel_level) - 1))
	{
	  if (pos [i] > 0)
	    res [i] = (pos [i] >> -rel_level) + 1;
	  else
	    res [i] = pos [i] >> -rel_level;
	}
      else
	res [i] = pos [i] >> -rel_level;
	
    return res;
  }

  template<typename aa_type>
  blitz::RectDomain<aa_type::Dim>
  aa_convert_to_level_coord (typename aa_type::cell_index& pos, int dest_level, int src_level)
  {
    typedef typename aa_type::IndexType IndexType;

    blitz::RectDomain<aa_type::Dim> coord;
    int rel_level = dest_level - src_level;
    
    if (!rel_level)
      {
	IndexType p ((pos.node () << dest_level) + pos.val ());

	return blitz::RectDomain<aa_type::Dim> (p, p + 1);
      }

    if (rel_level > 0)
      {
	IndexType p ((pos.node () << dest_level) + (pos.val () << rel_level));
	IndexType p1 (p + (1 << rel_level));
	
	return blitz::RectDomain<aa_type::Dim> (p, p1);
      }

    IndexType p ((pos.node () << dest_level) + (pos.val () >> (-rel_level)));
    return blitz::RectDomain<aa_type::Dim> (p, p + 1);
  }

  template<typename aa_type>
  blitz::RectDomain<aa_type::Dim>
  aa_convert_to_level_coord (aa_type* array, typename aa_type::cell_index& pos, int level)
  {
    typename aa_type::node_iterator it (array->at (pos.node (), false));

    return aa_convert_to_level_coord<aa_type> (pos, level, it.level ());
  }

  template<typename aa_type>
  typename aa_type::IndexType
  aa_min_rel_pos (typename aa_type::cell_index& p1, typename aa_type::cell_index& p2,
		  aa_type* array, int level)
  {
    typedef typename aa_type::IndexType IndexType;

    blitz::RectDomain<aa_type::Dim> d1 (aa_convert_to_level_coord<aa_type> (array, p1, level));
    blitz::RectDomain<aa_type::Dim> d2 (aa_convert_to_level_coord<aa_type> (array, p2, level));

    return rect_dom_min_rel_coord<aa_type::Dim> (d1, d2);
  }

  template<typename aa_type>
  typename aa_type::IndexType
  aa_min_rel_pos (typename aa_type::cell_index& p1, int level1, typename aa_type::cell_index& p2,
		  aa_type* array, int level)
  {
    typedef typename aa_type::IndexType IndexType;

    blitz::RectDomain<aa_type::Dim> d1 (aa_convert_to_level_coord<aa_type> (p1, level, level1));
    blitz::RectDomain<aa_type::Dim> d2 (aa_convert_to_level_coord<aa_type> (array, p2, level));

    return rect_dom_min_rel_coord<aa_type::Dim> (d1, d2);
  }

  template<typename aa_type>
  void 
  aa_sphere_neighbour (std::vector<typename aa_type::cell_index>& neighbours, 
		       aa_type* array, typename aa_type::cell_index& center, int r)
  {
    typedef typename aa_type::IndexType IndexType;

    typename aa_type::Node tmp_node (array->refine_level ());

    blitz::Array<bool, aa_type::Dim> n_map (IndexType (2*r + 1));
    blitz::RectDomain<aa_type::Dim> rad (-1, 1);
    
    std::deque<typename aa_type::cell_index> it_queue;

    n_map = false;
    n_map (r, r, r) = true;
    it_queue.push_back (typename aa_type::cell_index 
			(center.node (),
			 aa_convert_to_level_coord<aa_type::Dim> (center.val (), array->at (center.node (), false).level (), 
								  array->refine_level (), true)));

    while (!it_queue.empty ())
      {
	AAVirtNeighbourIterator<aa_type> it (array, it_queue.front (), &tmp_node, rad);
	it_queue.pop_front ();
	
	for (; !it.end_p (); ++ it)
	  {
	    typename aa_type::cell_index tmp_p (it.global_position (), it.local_position ());
	    IndexType rp (aa_min_rel_pos<aa_type> (center, array->refine_level (), tmp_p, array, array->refine_level ()));
	    IndexType p (rp + r);

	    if (rect_dom_in_p<aa_type::Dim> (n_map.domain (), p) && !n_map (p))
	      {
		n_map (p) = true;
		if (blitz::dot (rp, rp) <= r*r)
		  {
		    it_queue.push_back (tmp_p);
		    neighbours.push_back (tmp_p);
		  }
		
	      }
	  }
      }
    
  }

}; // namespace dbc

#endif
